%% OSL Task-postive ICA Example using
%% OHBA's ICA pipeLine (OIL)

% The OIL Pipeline tutorial is as follows:
% 1.) Beamform 4 subjects of CTF 2-back data
% 2.) Envelope estimation and down-sampling foolowed by spatial smoothing. 
% WE HAVE PROVIDED YOU WITH STAGES 1 & 2 AS THEY ARE SLOW TO COMPUTE.
% You can run this in your own time if you wish.
%------------------------------------------------------------------------

% YOU WILL CARRY OUT THE FOLLOWING:
% 3). Concatenation
% 4). Temporal ICA
% 5). 1st Level Stats
% 6). Group Stats
% 7.) Do Multiple Comparisons Corrections
% 8.) Generate Nifti Maps

% 9.) OPTIONAL extra: Rerun the ICA using ICASSO.

% You will need the CTF 2-back data
% and OSL 1.2

% Follow the script cell by cell.
% It should run without any changes so make to read the comments. 

% Henry Luckhoo
% henry.luckhoo@trinity.ox.ac.uk
% 13.08.12

% References:
% Brookes et al. 2011 PNAS - Resting State Analysis
% Luckhoo et al. 2012 Neuroimage - Task-postive analysis

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% Make sure that fieldtrip and spm are not in your matlab path
% SET THE BELOW LINE TO THE OSL DIRECTORY

osldir=['/home/hluckhoo/matlab/OSL/osl1.2.beta.4'];  % YOU NEED TO CHANGE THIS PATH
addpath(osldir);
osl_startup(osldir);
rmpath('/opt/analysis/spm8'); % If you have other versions of SPM on your path you should remove them too!

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

datadir='/home/hluckhoo/OSL_development_scripts/ICA_workshop/'; % directory where the data is - YOU NEED TO CHANGE THIS PATH

workingdir=[datadir '/spmfiles/epoched_data']; % this is the directory where the SPM files will be stored in
cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

spm_files={[workingdir '/espm8_ctf_2back_1'];
           [workingdir '/espm8_ctf_2back_2'];
           [workingdir '/espm8_ctf_2back_3'];
           [workingdir '/espm8_ctf_2back_4']};
       
% spm_files={[workingdir '/espm8_ctf_2back_1'];
%            [workingdir '/espm8_ctf_2back_2'];
%            [workingdir '/espm8_ctf_2back_3'];
%            [workingdir '/espm8_ctf_2back_4']};

structural_files = {'/home/hluckhoo/data/brookes_data/nback/sub_AV/AVMRI.nii';
                    '/home/hluckhoo/data/brookes_data/nback/sub_EH/ELHCrg.nii';
                    '/home/hluckhoo/data/brookes_data/nback/sub_JB/JBMRI.nii';
                    '/home/hluckhoo/data/brookes_data/nback/sub_JH/JHMRI.nii';};
                
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SETUP BEAMFORMER AND ICA OAT

% This is an example step up script for the whole analysis. You should look
% at the different parameters being set here and in "osl_setup_oil_ica.m".

oil = [];
oil.paradigm = 'task';

% Source Reconstruction Parameters
oil.source_recon = [];
oil.source_recon.dirname = [workingdir '/4-8Hz_espm8_ctf_2back_1_4subs_beta.4'];
oil.source_recon.D_epoched = spm_files;                                      % Here we are running all four subjects. 

oil.source_recon.mri = structural_files;
oil.source_recon.conditions = {'non-target' 'target'};
oil.source_recon.time_range=[-1 1];                                        % Time range of interest: here we have 2 second trials from -1s to 1s. If you set S2.time_range = '' the whole epoch is used. 
oil.source_recon.freq_range=[13 30];                                        % Frequency band of interest in Hz
oil.source_recon.modalities = {'MEGGRAD'};
oil.source_recon.fid_label.nasion='nas';
oil.source_recon.fid_label.lpa='lpa';
oil.source_recon.fid_label.rpa='rpa';
oil.source_recon.forward_meg ='MEG Local Spheres';

oil.source_recon.pca_dim = 0;
oil.source_recon.work_in_pca_subspace = 0;
oil.source_recon.force_pca_dim = 0;
oil.source_recon.gridstep = 6;

oil.enveloping = struct;
oil.enveloping.window_length = 0.5;

% Concatentation Parameters
oil.concat_subs = struct;
oil.concat_subs.normalise_subjects = 1; % Do we want to normalise subjects in general?

% ICA Parameters
oil.ica = struct;
oil.ica.num_ics = 20;   % We will do a 20 component decomposition today.

% First Level Parameters
oil.ica_first_level = struct;
% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to create the (num_regressors x num_trials) design matrix:
Xsummary={};Xsummary{1}=[1 0]; Xsummary{2}=[0 1]; 
oil.ica_first_level.design_matrix_summary=Xsummary;

% Contrasts to be calculated:
oil.ica_first_level.contrast={};
oil.ica_first_level.contrast{1}=[1 0]';       % 2-back non-targets
oil.ica_first_level.contrast{2}=[0 1]';       % 2-back targets
oil.ica_first_level.contrast{3}=[-1 1]';      % targets - non-targets 

% Group Level Parameters
oil.ica_group_level = struct;

% Check the OIL settings for any errors, missing values etc.
oil = osl_check_oil(oil);

%% 1-2.) RUN OAT stages 1 & 2

% If we were running the analysis for the first time through, we would
% uncomment and run the next two lines. However this could take some time
% even on fast computers. Therefore, the beamforming & source space
% enveloping has been done for you.wiki

oil.to_do = [1 1 1 1 1 1];        % Calls the beamformer "oil.to_do(1) = 1" and enveloping "oil.to_do(2) = 1".
oil = osl_run_oil(oil);


%%

% 
% 
% 
% S=[];
% S.paradigm='task';                                          % This is a critical flag to inform OSL whter to treat the data as "task" or "rest".
% S.oil_name = [workingdir '/4-8Hz_espm8_ctf_2back_1_4subs'];  % You would need to save a new OAT for each frequency band of interest
% S.D=spm_files;                                      % Here we are running all four subjects. 
% 
% % Source Reconstruction Parameters
% S.trigger={{'non-target'},{'target'}};
% S.freq_range=[13 30];                                        % Frequency band of interest in Hz
% S.time_range=[-1 1];                                        % Time range of interest: here we have 2 second trials from -1s to 1s. If you set S2.time_range = '' the whole epoch is used. 
% S.recon_method='beamform';
% S.gridstep=6;                                               % in mm, using a lower resolution here than you would normally, for computational speed
% S.modalities={'MEGGRAD'};                                 % Axial gradiometers in use
% S.fsample=600;                                              % CTF data has a sampling frequency of 600Hz
% S.pca_dim = 0;                                              % Uses full rank covariance matrix                                                            
% S.mri={'/home/hluckhoo/data/brookes_data/nback/sub_AV/AVMRI.nii';
%        '/home/hluckhoo/data/brookes_data/nback/sub_EH/ELHCrg.nii';
%        '/home/hluckhoo/data/brookes_data/nback/sub_JB/JBMRI.nii';
%        '/home/hluckhoo/data/brookes_data/nback/sub_JH/JHMRI.nii';};
% S.fid_label.nasion='nas';
% S.fid_label.lpa='lpa';
% S.fid_label.rpa='rpa';
% % ICA Parameters
% S.normalise_subjects=1;                                     % Do we want to normalise subjects in general?
% S.normalise_vox=0;                                          % Do we want to normalise voxels before the ICA?
% S.hilbert_downsampling=0.5;                                 % This is the window length for the down-sampling. 0.5s gives 4 timepoints per 2s trial.
% S.num_ics=20;                                               % This is the number of ICs sought
% 
% % Stats Parameters
% % Xsummary is a parsimonious description of the design matrix.
% % It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% % is a condition no. This will be used (by expanding the conditions over
% % trials) to create the (num_regressors x num_trials) design matrix:
% Xsummary={};Xsummary{1}=[1 0]; Xsummary{2}=[0 1]; 
% S.design_matrix_summary=Xsummary;
% 
% % Contrasts to be calculated:
% S.contrast={};
% S.contrast{1}=[1 0]';       % 2-back non-targets
% S.contrast{2}=[0 1]';       % 2-back targets
% S.contrast{3}=[-1 1]';      % targets-non-targets 
% 
% % Build the OAT
% oil = osl_setup_oil(S);


%% 3-6.) Run stages 3. concatenation of subjects, 4.) ICA decomposition, 5.) first level stats & 6.) group level stats.

% Here you load in an existing OAT result.
clear oil
oil = osl_load_oil([workingdir '/4-8Hz_espm8_ctf_2back_1_4subs' '.oil']);        % Here we load in a pre-existing OAT result. 
oil.to_do = [0 0 1 1 1 1];                                                                                 

oil = osl_run_oil(oil);

%%%%%%%%%%%%%%%%%%%
%% 7. Significance of the ICs
clear oil
oil = osl_load_oil([workingdir '/4-8Hz_espm8_ctf_2back_1_4subs' '.oil']);     % Load in the oil we have just run

con=3;                                                                      % Select the targets-non-targets contrast - this is our differential contrast. We also calculated the main effects (con = 1 & 2).
[~, order]=sort(abs(oil.ica_group_level.results.tstats(:,con)),'descend');  % We want to rank our components in order of descending significance

% Let's try weak multiple comparisons correction first.

pvals_corrected_ordered=correct_ica_pvals(oil.ica_group_level.results.pvals(order,con),0.05,'weak'); % Carry out multi-step-up correction (Nichols & Hayasaka, 2003)
disp('Weakly corrected p-values');
disp(pvals_corrected_ordered);

% Now let's try strong multiple comparisons correction.

pvals_corrected_ordered=correct_ica_pvals(oil.ica_group_level.results.pvals(order,con),0.05,'strong'); % Carry out Bonferroni correction.
disp('Strongly corrected p-values');
disp(pvals_corrected_ordered);

% Note how fewer of the components survive the strong correction. In this
% case the p-values are all very small so we can afford to apply strong
% correction. However in some cases strong correct will be too severe. In
% this case, a weaker control of the family-wise error rate can be adopted
% by using the multi-step-up test.

%%%%%%%%%%%%%%%%%%%
%% 8. Generate IC correlation maps Nifti files

[maps] = osl_save_nii_ica_maps(oil,'correlation',order);                         % Generate Nifti Correlation maps in order of significance
fslview(maps)

% On the fslview GUI:  set the positive AND negative colormaps using the "i"
% button and THEN set min to 0.3 and max to 0.7 and. Have a look at the various components. 
% NOTE that this analysis is only run with 4 subjects so the results won't look
% great.

% From this data we have previously found significant components in the
% frontal areas, the visual cortex. motor cortices and the SPL. Can you
% find any components that might match those regions? Are they significant?

[maps] = osl_save_nii_ica_maps(oil,'scaled_covariance',order);                         % Generate Nifti Correlation maps in order of significance
fslview(maps)
%%%%%%%%%%%%%%%%%%%
%% 9.) OPTIONAL: Rerun the ICA but using ICASSO instead.

oil=osl_load_oil([workingdir '/4-8Hz_espm8_ctf_2back_1_4subs' '.oil']); % Load in the oil we have just run
oil.ica.icasso_its=5;                                                       % NOTE THAT WE ARE ONLY RUNNING 5 ITERATIONS. This is too small. 30 would be a more sensible number for this analysis
                                                                            % We could edit any of the  parameters like this. e.g. num_ics or hilbert_downsampling .
oil.to_do = [0 0 0 1 1 1];
oil = osl_run_oil(oil);                                                     % Note that ICASSO is running multiple iterations - as a result it takes longer.


%% AGAIN LOAD THE RESULTS
clear oil
oil = osl_load_oil([workingdir '/4-8Hz_espm8_ctf_2back_1_4subs' '.oil']);     % Load in the oil we have just run

con=3;                                                                      % Select the targets-non-targets contrast - this is our differential contrast. We also calculated the main effects (con = 1 & 2).
[~, order]=sort(abs(oil.ica_group_level.results.tstats(:,con)),'descend');                                     % We want to rank our components in order of descending significance

% Let's try weak multiple comparisons correction first.

pvals_corrected_ordered=correct_ica_pvals(oil.ica_group_level.results.pvals(order,con),0.05,'weak'); % Carry out multi-step-up correction (Nichols & Hayasaka, 2003)
disp('Weakly corrected p-values');
disp(pvals_corrected_ordered);

% Now let's try strong multiple comparisons correction.

pvals_corrected_ordered=correct_ica_pvals(oil.ica_group_level.results.pvals(order,con),0.05,'strong'); % Carry out Bonferroni correction.
disp('Strongly corrected p-values');
disp(pvals_corrected_ordered);

% Note how fewer of the components survive the strong correction. In this
% case the p-values are all very small so we can afford to apply strong
% correction. However in some cases strong correct will be too severe. In
% this case, a weaker control of the family-wise error rate can be adopted
% by using the multi-step-up test.

%% GENERATE THE CORRELATION MAPS

[maps] = osl_save_nii_ica_maps(oil,'correlation',order);                         % Generate Nifti Correlation maps in order of significance
fslview(maps)

% On the fslview GUI:  set the positive AND negative colormaps using the "i"
% button and THEN set min to 0.3 and max to 0.7 and. Have a look at the various components. 
% NOTE that this analysis is only run with 4 subjects so the results won't look
% great.

% Do you think the ICASSO results look better? Again we have limited
% ourselves to only 4 subjects and only 5 iterations so we can't really
% tell with this data.



